Experimental!¶
Disclaimer: All pixel size extraction code is a first pass guess. I do not have access to a Varian Linac.
In [1]:
import numpy as np
import matplotlib.pyplot as plt
import IPython.display
import pydicom
In [2]:
!pip install pymedphys==0.19.0dev0
In [3]:
import pymedphys
import pymedphys._wlutz.core
import pymedphys._wlutz.reporting
import pymedphys._vendor.pylinac.winstonlutz
In [4]:
pymedphys.__version__
Out[4]:
In [5]:
def display_markdown(string):
IPython.display.display(IPython.display.Markdown(string))
display_markdown('## Example heading')
In [6]:
dicom_wlutz_paths = pymedphys.zip_data_paths('denis_wlutz_images.zip')
In [7]:
dicom_wlutz_paths
Out[7]:
In [8]:
dicom_files = [pydicom.read_file(str(path), force=True) for path in dicom_wlutz_paths]
In [9]:
dcm_header = dicom_files[0]
dcm_header
Out[9]:
In [10]:
dcm_header.ImagePlanePixelSpacing
Out[10]:
In [11]:
sad = float(dcm_header.RadiationMachineSAD)
panel_adjustment = -float(dcm_header.XRayImageReceptorTranslation[2])
panel_ssd = panel_adjustment + sad
guess_at_pixel_spacing_at_iso = np.array(dcm_header.ImagePlanePixelSpacing).astype(float) / panel_ssd * sad
dx, dy = guess_at_pixel_spacing_at_iso
dx, dy
Out[11]:
In [12]:
half_range_x = dcm_header.Columns * dy / 2
half_range_y = dcm_header.Rows * dx / 2
In [13]:
x = np.linspace(-half_range_x, half_range_x, dcm_header.Columns)
y = np.linspace(-half_range_y, half_range_y, dcm_header.Rows)
In [14]:
jaw_pos = {
coll.RTBeamLimitingDeviceType: np.array(coll.LeafJawPositions).astype(float)
for coll in dcm_header.ExposureSequence[0].BeamLimitingDeviceSequence
}
jaw_pos
Out[14]:
In [15]:
field_size_x = np.diff(jaw_pos['ASYMX'])[0]
field_size_y = np.diff(jaw_pos['ASYMY'])[0]
edge_lengths = [field_size_x, field_size_y]
In [16]:
penumbra = 2
bb_diameter = 8
In [17]:
gantries = np.array([
np.round(dcm.GantryAngle, 2) for dcm in dicom_files
])
gantries
Out[17]:
In [18]:
colls = np.array([
np.round(dcm.BeamLimitingDeviceAngle, 2) for dcm in dicom_files
])
colls
Out[18]:
In [19]:
images = [dcm.pixel_array for dcm in dicom_files]
scaled_images = [
img[::-1, :] / 2 ** 16 for img in images
]
In [20]:
for img in scaled_images:
plt.figure()
plt.imshow(scaled_images[0])
plt.colorbar()
In [21]:
def get_initial_centre(x, y, img):
wl_image = pymedphys._vendor.pylinac.winstonlutz.WLImageOld(img)
min_x = np.min(x)
dx = x[1] - x[0]
min_y = np.min(y)
dy = y[1] - y[0]
field_centre = [
wl_image.field_cax.x * dx + min_x,
wl_image.field_cax.y * dy + min_y
]
return field_centre
In [25]:
def analysis(x, y, img, collimator_angle):
field = pymedphys._wlutz.imginterp.create_interpolated_field(x, y, img)
initial_centre = get_initial_centre(x, y, img)
pymedphys._wlutz.reporting.image_analysis_figure(
x,
y,
img,
None,
initial_centre,
collimator_angle,
bb_diameter,
edge_lengths,
penumbra,
)
plt.title('Initial Centre')
field_centre, field_rotation = pymedphys._wlutz.findfield.field_centre_and_rotation_refining(
field, edge_lengths, penumbra, initial_centre, pylinac_tol=np.inf, fixed_rotation=collimator_angle
)
bb_centre = pymedphys._wlutz.findbb.optimise_bb_centre(
field, bb_diameter, edge_lengths, penumbra, field_centre, field_rotation, ignore_pylinac=True
)
pymedphys._wlutz.reporting.image_analysis_figure(
x,
y,
img,
bb_centre,
field_centre,
collimator_angle,
bb_diameter,
edge_lengths,
penumbra,
)
plt.title('PyMedPhys Basinhopping Method')
try:
pylinac = pymedphys._wlutz.pylinac.run_wlutz(
field, edge_lengths, penumbra, field_centre, collimator_angle)
pymedphys._wlutz.reporting.image_analysis_figure(
x,
y,
img,
pylinac['v2.2.6']['bb_centre'],
pylinac['v2.2.6']['field_centre'],
collimator_angle,
bb_diameter,
edge_lengths,
penumbra,
)
plt.title('Pylinac v2.2.6 Filter and Profile Method')
pymedphys._wlutz.reporting.image_analysis_figure(
x,
y,
img,
pylinac['v2.2.7']['bb_centre'],
pylinac['v2.2.7']['field_centre'],
collimator_angle,
bb_diameter,
edge_lengths,
penumbra,
)
plt.title('Pylinac v2.2.7 Filter and Scikit-Image Method')
except Exception as e:
print(e)
return np.array(field_centre) - np.array(bb_centre)
In [30]:
for img, coll, gantry in zip(scaled_images, colls, gantries):
display_markdown(f'## Gantry {coll} | Collimator {gantry}')
deviation = np.round(analysis(x, y, img, coll), 2)
display_markdown(
f'PyMedPhys field centre - BB centre (mm):\n\n```python\n[x, y] = [{deviation[0]}, {deviation[1]}]\n```'
)
plt.show()
